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Abstract 


An  overview  is  presented  of  the  authors'  recent  theoretical  and 
experimental  results  on  reliable  and  computationally  efficient  a-posteriori 
error  bounds  for  finite  elanent  solutions.  These  estimates  are  composed 
from  error  indicators  evaluated  on  the  individiial  elements,  and  these 
indicators  in  turn  allow  for  a very  effective  approach  to  the  effective 
construction  of  optimal  meshes.  Finally,  some  views  are  presented  about 
possible  future  trends  in  the  development  of  finite  element  software  and 
an  outline  is  given  of  the  design  of  an  experimental  finite  element  system 
currently  under  development  which  incorporates  many  of  these  ideas  and 


results. 


1 . Introduction 


In  recent  years  considerable  progress  has  been  made  in  the  theoretical 
analysis  and  the  variety  and  depth  of  applications  of  the  finite  element 
method  (see,  e.g.,  [1  ]).  The  pace  of  these  developments  does  not  appear 
to  be  diminishing  (see,  e.g.,  [2]).  It  is  surprising,  however,  that  a 
subject  of  considerable  practical  importance,  namely,  the  assessment  of  the 
reliability  of  the  results  of  the  finite  element  computations,  has  not 
been  studied  more  intensively. 

Many  of  the  current  reliability  studies  consist  of  various  comparisons 
for  different  benchmark  problems.  Theoretical  error  estimations  are  almost 
always  of  an  a-priori  type  and,  as  inqxsrtant  as  such  estimates  are  for  the 
theory  of  the  method,  they  are  rarely  of  much  use  for  practical  error 
determinations.  In  particular,  since  in  practical  applications  an  accuracy 
of  5-10%  may  be  entirely  acceptable,  any  error  bound  that  typically  repre- 
sents an  overestimate  by  a factor  of  five  is  clearly  worthless. 

It  appears  that  for  any  effective  computation  of  reliable  error  bounds 
for  finite  element  solutions  we  should  use  a-posteriori  estimates  that  are 
based  on  information  obtained  during  the  solution  process  itself.  Moreover, 
such  estimates  may  also  provide  a tool  for  the  adaptive  design  of  optimal 
finite  element  meshes. 

This  is  the  topic  of  this  paper.  More  specifically,  we  present  here 
an  overview  of  recent  theoretical  and  experimental  results  ([3] , [4] , [5] , [6]) 
on  reliable  and  conputationally  efficient  a-posteriori  error  bounds  for  the 
finite  element  method.  These  estimates  are  composed  from  "error- indicators" 
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evaluated  on  the  elanents  used  in  the  solution,  and  these  error  indicators 
in  turn  allow  for  a very  effective  approach  to  the  adaptive  construction 
of  optimal  meshes.  After  a brief  discussion  in  Section  2 of  typical  a- 
priori  estimates  and  their  shortcomings  for  practical  purposes,  we  present 
in  Section  3 some  of  the  principal  results  about  the  a-posteriori  estimates. 
This  includes  also  a few  illustrative  numerical  experiments.  Then  Section  4 
addresses  the  design  of  optimal  finite  element  meshes  using  these  estimates. 
Finally,  in  Section  5 we  sketch  some  views  about  possible  future  trends  in 
the  development  of  finite  element  software  and  outline  an  experimental 
finite  element  system  currently  under  development  which  incorporates  many 
of  these  ideas  and  results. 
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2,  Rates  of  Convergence 

In  this  section  we  present  some  typical  exan^les  of  theorems  about  the 

rate  of  convergence  of  the  finite  element  solutions  to  the  exact  solution. 

For  reasons  of  sitrplicity  we  restrict  the  discussion  to  the  sijnple  model 
2 

problem  in  R 


(2.1) 


■Au  + u = f,  on  c R'^ 


= 0,  on  952 


where  52  is  an  L- shaped  domain  as  shown  Figure  1 

in  Figure  1,  n represents  the  outward 

normal  of  352,  and  f is  a suitable  given  function  on 

On  S we  consider  a family  P of  quasi-uniform  triangular  meshes  A; 
that  is,  the  ratio  of  the  largest  triangle-side  to  the  smallest  triangle- 
side  occurring  in  A is  bounded  by  a constant  depending  only  on  P.  For 
any  A € P let  S(A)  be  the  set  of  all  continuous  functions  on  which  are 
linear  on  each  triangle  of  A.  The  dimension  of  the  space  S(A)  shall  be 
denoted  by  N(A).  In  addition,  let  Qp  be  the  set  of  all  polynomials  of 
degree  at  most  p on  2 and  Np  the  dimension  of  Qp. 

For  any  A € P the  finite  element  solution  of  (2.1)  is  now  the  function 


u = u(A)  € S(A)  such  that  for  all  v € S(A) 
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(2.2) 


1 

2 


' 8U  dv  ^ au  dV 
axj^  3X2  3X2 


+ 


uv  dx  = 


/ fvdx  . 


Correspondingly  u = u(p)  e is  the  function  for  which  (2.2)  holds  for  all 
V € Qp.  The  exact  solution  of  (2.1)  is  denoted  by  Uq,  and  we  measure  the 
errors  e(A)  = Ug  - u(a)  and  e(p)  = Ug  - u(p)  in  the  energy  norm 


(2.3) 


- / 


+ e 


dx 


The  rate  of  convergence  depends  on  the  smoothness  of  the  exact  solution 
Ug.  Mathematically,  it  is  advantageous  to  introduce  the  so-called  Besov 
space  B?  (ffi)  (see,  e.g.,  [ 7])  and  for  u«  € b“  (ffi)  to  express  the  error 

U A. 

in  terms  of  the  norm  of  that  space.  We  shall  not  enter  into  details  here, 

but  note  that  H°(S2)  c B?  c H°^^(2)  where  H^(S2)  is  the  usual  Sobolev 

space  of  fractional  order  a on  S2.  We  mention  also  that  the  function 
2/3 

1 cos  20/3- -representing  the  natural  singularity  of  our  problem  on  the  L- 
shaped  domain  S2--belongs  to  b5^^(S2)  but  not  to  H^'^^(ffi). 

For  the  following  rate  of  convergence  result  we  consider  sequences  of 
meshes  € P for  which  S(A^)  c and  N(Aj)  -♦  «,  as  i -*■  <».  We  call 

this  a sequence  of  proper  refinements. 

Theorem  2.1;  Stqppose  that  the  exact  solution  Ug  of  (2.1)  belongs  to 

02  «(2),  then 

1-a 

||e(p)||g  < C llugll  ^ , as  Np  ^ 

®2,- 


(2.4) 


OP 


i 
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where  the  constant  C is  independent  of  p and  but  depends  on  a 
and  , etc.  If  1 < a < 2,  then  for  any  sequence  of  proper  refinements 

1-a 

(2.5)  ||e(A)||E  2 C N(a)~^  liupll  ^ , as  N(a)  - », 

Bo 

Z 

^ 9 

wheie  C is  independent  of  Uq  and  the  meshes  used,  but  once  again  depends 
on  a,  2,  etc. 

For  the  finite  element  solution  the  case  a = 2 requires  special  atten- 
tion. But  we  note  that  when  linear  finite  elements  are  used,  the  rate  cannot 
-1/2 

be  better  than  N ' . 

It  is  important  that  conversely  the  rate  of  convergence  completely 
characterizes  the  regularity  of  the  solution  Uq: 

Theorem  2.2:  Suppose  that  either 

1-a 

(2.6)  l|e(p)||p  < K , as  N , 

c p p 

or  for  1 < a < 3/2  and  for  any  sequence  of  proper  refinements  , 

1-a 

(2.7)  l|e(A)||g  ^ K NU)~^  , as  N(A)  ^ 

then  Uq  € B2  ^(2)  and 

(2.8)  IIuqII  ^ s C[K  + (/  u^dx)^/2]  . 

®2,-  ^ 

In  particular,  since  r^^^  cos  26/3  € b5^^(2),  the  natural  singularity 

slows  down  the  convergence  rate  to  [*]  as  compared  to  the  value 

■*  -1/3+e 

Usually  only  the  value  N ' is  proved. 


^-1/2 
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for  smooth  solutions.  In  [ 8 ]--and,  with  more  details  in  [ 9 ] 

--it  was  shown  that  for  properly  refined  sequences  of  meshes  the  rate 
-1/2 

N can  be  obtained  even  in  the  presence  of  the  natural  singularity. 

In  tliat  case  it  is  natural  to  use  weighted  Besov  spaces  with  weights  that 
depend  on  the  refinement.  This  type  of  result  reflects  the  importance  of 
the  use  of  properly  refined  meshes. 

The  conparison  between  the  behavior  of  the  polynomial  approximations 
Up  and  the  finite  element  solutions  u(a)  shows  that  for  the  same 
degree  of  freedom  N = Np  = N(a)  the  polynomials  give  the  same  rate  of  con- 
vergence if  the  solution  is  badly  behaved.  Of  course,  for  smooth  Uq 
the  convergence  rate  for  the  polynomial  solutions  will  be  better  or  at  least 
never  worse  than  that  of  the  u(A).  This  is  the  basis  of  the  p-convergence 
approach  developed  in  [10],  [11]. 

For  practical  purposes  the  estimates  of  Theorem  2.1  are  not  useful. 

The  constants  are  computationally  unavailable,  and  so  is  the  norm  of  the 
exact  solution.  The  estimates  (2.4),  (2.5)  are  asymptotic  in  nature  and 
hence  characterize  the  convergence  rate  for  large  N.  For  smaller  N,  the 
decrease  of  the  error-norm  with  growing  dimension  is  often  faster  than  the 
asynqjtotic  rate.  This  effect  has  been  observed  in  many  practical  situations. 
It  is  intuitively  understandable  since  changes  of  Uq  have  to  be  compared  to 
the  mesh  size. 

As  stated  at  the  outset,  the  discussion  in  this  section  is  only  meant  to 
illustrate  a typical  situation.  Theorems  of  the  form  given  here  apply  also  in 
many  more  general  cases  (see,  e.g.,  [12]). 
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3.  A-Posteriori  Error  Estimates 

The  discussion  of  the  previous  section  shows  that  a-priori  estimates, 
such  as  (2.4), (2.5),  are  not  directly  usable  for  practical  error  assess- 
ments. In  order  to  obtain  conputable  and  reliable  bounds,  we  need  to  go 
to  a-posteriori  estimates  which  use  information  obtained  during  the  solu- 
tion process  itself.  The  theory  of  these  estimates  differs  somewhat  for 
1 d 

problems  in  R and  those  in  R with  d > 1,  and  accordingly  we  discuss 
these  cases  separately. 

3.1  One -Dimensional  Problems 

Once  again- -for  sin5)licity--we  restrict  the  presentation  to  a sin^jle 
model  problem.  The  results  also  extend  to  more  general  cases,  as,  for 
instance,  the  beam  equation.  Consider  the  two-point  boundary  problem 

' ^ ^ * b(x)u  = f(x),  0 < X < 1 

(3.1) 

u(0)  = u(l)  = 0 


with  sufficiently  smooth  coefficients  a,b,f  on  [0,1]  such  that  a(x)  > a > 0, 
and  b(x)  > 0 for'O  5 x £ 1. 

Let  P be  a family  of  partitions 

(3.2)  a:  0 = Xq  < x^  < . . . < x^_j  < x^  = 1,  m = m(A)  > 1 , 

on  the  interval  [0,1],  and  introduce  the  notations 
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h.(A)  = ,,  j = 1, 

J 3 J-1 


. . . .m 


(3.3) 


h(A)  = max  h.(A),  h(A)  = min  h.(A)  . 
j=l m ^ j=l,...,m  ^ 


The  family  P will  be  called  ic-regular  if  there  are  constants  X > 0,  <5 
such  that 


(3.4) 


h(A)  > Xh(A)  , V A € P 


For  this  presentation  we  shall  use  linear  elements  on  the  intervals 
I^(A)  = [xjjjXj],  j = l,...,m,  to  solve  (3.1).  As  before,  we  denote  by 
u(A)  the  resulting  finite  element  solution  and  by  Uq  the  true  solution 
of  our  problem.  On  every  separate  interval  1^  (A)  the  function  u(A)  is 
linear  and  hence  the  residuals 


(3.5) 


rj(x)  = 


a(x)  u(A)(x)  + b(x)u(A)  (x)-f(x) , Xj_j^  < x s x ^ 


J 1 , . . . ,m, 


and  their  integrals 


(3.6) 


X. 

.2  _ 


Vj(A)  = / rj(x)'^dx,  j = 1 m 

""j-l 


are  easily  computable.  We  introduce  the  error  indicators 

,2.  1 


(3.7) 


“ T2  — r 

ac|  Cxj.j.x.)) 


9 J 1 > • • • >in 


and  the  error  estimator 


(3.8) 


- 9 - 

e(A)^  = I e.(A)^  . 
j=l  J 

Then  the  following  result  holds;  see  [5] , [6] : 

Theorem  3.1:  Let  Uq  be  the  true  solution  of  (3.1)  and  u(A)  the 
finite  element  solution  for  any  A € P.  Then  without  any  regularity  assump- 
tion about  P the  error  e(A)  = Uq  - u(A)  satisfies 

(3.9)  ||e(A)||p^  S ^ 6(A)(l+0(h(A))),  as  h(A)  0 . 

The  bound  of  the  0-term  depends  on  a,b  but  not  on  Uq.  If  P is  K-regular 
with  1 i K < 2 and  Uq(x)  has  only  simple  roots  in  [0,1],  then 

(3.10)  l|e(A)||p  = 6(A)(l+0(h(A)'")),  as  h(A)  - 0 , 
where  k = 1 - k/2. 

In  the  present  case  the  energy  norm  used  in  (3.9)  and  (3.10),  of 
course,  has  the  form 

(3.11)  \\e\\l  = / (a(x)(^)^  + b(x)e2)dx  . 

Note  that  in  the  case  of  (3.10)  the  error  indicator  (3.8)  provides  the 
exact  error  for  h -*•  0.  In  other  words,  the  effectivity  con5tant 

l|e(A)|L 

(3.12) 

tends  to  one  for  h 0.  Even  the  general  estimate  (3.9)  ensures  that 
asymptotically  we  have  0(A)  s VTZ/^  < 1.103. 
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I 


In  applications  we  may  readily  accept,  say,  1/2  5 0(A)  5 2,  that  is 
a 1001  error  in  0(A)  with  respect  to  the  limiting  value  0(A)  =1.  On 
the  other  hand,  a corresponding  error  in  the  solution  relative  to  u^,  that 
is,  1/2  s ||e(A) llg/llugllj;  s 2,  would  be  coirpletely  unacceptable.  This  points 
out  once  again  the  essential  difference  in  the  estimates  in  this  section  and 
those  of  Section  2. 

As  an  illustration  of  the  effectivity  of  Theorem  3.1,  we  consider  the 
following  example: 

' ^ ^ (x+a)*^!  = f,  0 < X < 1,  a > 0, 

(3.12) 

u(0)  = u(l)  = 0 . 

Here  f was  chosen  such  that  the  exact  solution  of  (3.12)  is 


(3.13)  Uq(x)  = (x+a)^  - [a^(l-x)+(l+a)^x]  . 


Note  that  for  small  a and  negative  r severe  near  singularities  may  be 
created.  In  the  numerical  examples  given  below  we  chose 

(3.14)  P = 0,  q=l,  » 

in  which  case  IIuqH^  = 6.09811.  Two  different  types  of  meshes  are  used, 
namely,  uniform  meshes  with  h^  = 1/m  and  asyiUJtotically  optimal  meshes  in 
the  sense  of  Section  4 below.  Table  1 below  shows  the  coordinates  of  the 
nodal  points  of  the  asymptotically  optimal  mesh  with  m = 10. 


• f V • 
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Table  1 


j 

Asymptotically  Optimal  Mesh  for  m = 

10 

X. 

I 

X. 

1 

j 

j 

0 

0 

4 

.01443 

8 

.11781 

1 

.00207 

5 

.02398 

9 

.26831 

2 

.00487 

6 

.03732 

10 

1 

3 

.00877 

7 

. 06318 

In  Table  2 we  show  the  relative  error  E = 100  ||e||g/||uQ||g  in  percent 
of  the  energy  norm  of  the  exact  solution  and  the  effectivity  constant  9 
of  (3.12)  for  the  two  types  of  meshes.  In  addition,  two  quantities 
CO  and  Eq  are  included  which  relate  to  the  optimality  of  the  mesh  as 
discussed  in  Section  4 below.  More  specifically,  co  is  the  ratio  of  the 
maximal  and  minimal  error  indicators.  As  we  shall  see  in  Section  4,  a 
mesh  is  asymptotically  optimal  if  all  its  error  indicators  are  equal. 

Thus  CO  is  a measure  of  optimality.  The  value  Eq  is  the  value  of  E 
for  the  asymptotically  optimal  partition. 

The  reliability  and  effectivity  of  the  estimates  are  evident.  A 26% 
error  relative  to  the  norm  of  the  exact  solution  corresponds  only  to  a 
15%  error  in  estimating  the  actual  error  by  means  of  the  estimator  (3.12). 
For  the  "more  optimal"  meshes  the  results  are  even  better. 


O-..- 
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Table  2 


Uniform  mesh 


m 

E 

Eo 

0 

CO 

5 

85.301 

22.613 

.1706 

8. 84  (+6) 

10 

73.768 

11.306 

.2950 

2. 34 (+7) 

20 

58.784 

5.653 

.4702 

5. 31 (+7) 

40 

41.933 

2.827 

.6708 

1.11 (+8) 

80 

26.316 

1.413 

.8419 

2. 19 (+8) 

Asymptotically  optimal  mesh 

m E 

Eq 

e 

CO 

5 

22.243 

22.613 

.6524 

5.854 

10 

11.289 

11.306 

.9025 

2.274 

20 

5.652 

5.653 

.9757 

1.372 

40 

2.826 

2.827 

.9940 

1.111 

80 

1.413 

1.413 

.9984 

1.031 

\ 

I 


n 
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3. 2 Higher- Dimensional  Problems 

Once  more  we  restrict  ourselves  to  a special  model  problem,  namely, 
the  two-dimensional  boundary  value  problem 
2 

■ ^ ^ ^ + bu  = f(x),  V X € S2 

i,j=l  ^''i 


(3.15) 


uCx)  = g(x),  V X € aS2  . 


Here  the  coefficients  a^j,b  are  constants  with 


(3.16) 


®ij^i^j  - i(x^+X2),  V X € 3^2  " ^21 


b > 0 , 


f is  continuous  on  ffi  and  continuously  differentiable  on  S2,  and  g is 

continuous  and  piecewise  continuously  differentiable  on  aS2. 

2 

IVe  restrict  the  domain  S2  to  sets  of  R bounded  by  lines  that  are 
parallel  to  the  coordinate  axis.  For  example,  ' 

Q may  be  the  L- shaped  domain  of  Figure  1 in  > m oq  i m i 

Section  2. 

I h»-i  1-  ■ 

The  admissible  meshes  A on  shall 

il — <—ii 

be  sets  A = {qj } of  squares  qj  of  the  forni  , , ifa  | j 


(3.17) 


such  that 


qj  = (x  € R'^  a < X^  < a + h,  i = 1,2},  h = h(qj)  > 0 , 
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(i)  for  any  qj  € A the  closure  belongs  to 


(3.18) 


(ii)  U q.  = a; 

q^€A  ^ 


(iii)  q^  n qj  = 0 for  any  q^,qj  € A,  i j ; that  is,  q^  n q^ 
consists  only  of  boundary  points  of  q^^  and  q^ . 

A typical  mesh  on  a square  ffi  is  given  in  Figure  2.  The  comer  points  of 
the  squares  q^.  € A are  the  nodal  points  of  the  mesh.  A nodal  point  v d S 
is  regular  if  it  is  a comerpoint  of  all  the  squares  qj  € A for  which 
V € q j ; otherwise  v is  irregular.  Note  that  then  any  nodal  point  on 
aa  is  necessarily  regular.  In  Figure  1 the  regular  and  irregular  nodes 
are  marked  by  black  dots  and  open  squares,  respectively. 

A collection  P of  admissible  meshes  A on  a is  N-regular  if  for 
any  A € P the  number  of  irregular  nodes  on  any  side  of  a square  qj  d A 
is  at  most  eqioal  to  W. 

For  a given  admissible  mesh  A on  a let  5(A)  be  the  set  of  all 
continuous  functions  u = u(xj,X2)  on  S which  are  bilinear  in  Xj^,X2  on 
every  q^  d A.  Obviously,  any  u € 5 (A)  is  then  uniquely  determined  by  its 

O 

values  in  the  regular  nodes  of  A.  By  5 (A)  we  denote  the  subset  of  all 
functions  of  5 (A)  which  are  equal  to  zero  on  aa. 

In  order  to  facilitate  the  error  estimation  on  the  boundary  we  add  the 
technical  condition  for  g that  for  any  qj  € A v4iich  has  a side  a on  aa 
we  have  with  a suitable  constant  c > 0 


•o-'J-.V 
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(3.19) 


/ (g'-u(i)')^dx  < / (g-u(A))^dx,  a = aq.  fl  dS2 

a CT  ■’ 


This  holds,  for  instance,  for  quadratic  polynomials  g. 

For  a given  admissible  mesh,  the  finite  element  solution  of  (3.15)  is 
now  the  (unique)  function  u(a)  € S(a)  for  which 


(3.20) 


/ ( I ^ ^ V V € S(A) 

2 i,^l  ^^i  2 


and  u(A)  = g holds  at  all  nodal  points  of  A on  32.  As  before  Uq  denotes 
the  exact  solution  of  (3.15)  and  we  are  interested  in  estimating  the  error 
e = Uq  - u(a)  in  the  energ>'  norm 

As  in  the  one -dimensional  case  we  may  conpute  on  each  q^  € a the  residual 


(3.21) 


rj(x) 


l — a..  — 

1,5=1  ”1  “i 


U(A)  + bu(A)  - f(x),  V X € q. 


If,  for  instance,  = a22  = 1»  Si-^^2  = 0,  b = 0 then  it  is  easily  seen 

that  r.(x)  = 0 in  any  q.  for  which  f(x)  = 0,  x € q. . This  prevents  us 
J 1 j 

from  proceeding  analogously  to  the  one -dimensional  case.  Instead  we  use  an 
approach  which  is  equivalent  to  that  presented  in  [ 3 ] . 

Let  qj  € A be  any  square  of  the  mesh  

and  as  indicated  in  Figure  3 let  0^,02 .^2 

be  the  four  sides  of  qj . For  a square  qj  * ' 

in  the  interior  of  the  domain,  we  denote  by 
J and  J-  the  junps  of  u(A)  on  a-i  e: 

and  Oj^,  respectively,  and  correspondingly.  Figure  3 


- 16  - 


by  J and  J-  those  of  — - u(a)  on  Oo  and  a,,  respectively.  Now 
we  con5)ute  the  following  quantities: 


If  c 2 then 


r-.  r .2  1 h(qj' 

(l)  TiCq^)  - ^ ^ 


/ r- (x)^dx 

n J 


(3.22) 


(ii)  P^C^j)  = 


^ll^C^i)  / J^(x)^dx,  a = a^,a^ 

fl 


1 2 
IT  a22h(qi)  / J^(x)  dx,  a = 02, c 
^ a 


and  the  error  indicator  of  q^  is  defined  by 

(3.23)  e(q.)2  = q(q.)2  + p^(q.)2  ^ p^(q.^2  ^ ^ 

On  the  other  hand,  if  q^  has  sides  in  common  with  32  then  for  any 
such  side  a € q^  n 32  the  corresponding  p-term  in  (3.23)  is  replaced  by 

(3.24)  pCT^^^j)  “ TTi^y  ^ (g‘“CA))^<lx  • 

Here  and  in  (3.22)(i),  a and  a are  the  constants  of  (3.16). 

As  in  the  one- dimensional  case  we  introduce  now  the  error  estimator 


(3.25) 


(A)  « I e(q,r  . 


Then  the  following  result  holds  [ 3 ] : 
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Theorem  3.2:  The  error  e * Uq  - u(A)  between  the  exact  solution  and 
the  finite  element  solution  of  (3.15)  satisfies 

(3.26)  i e(A)  < llellg  < K e(A) 

with  some  K > 1 which  is  independent  of  A but  depends  on  a,  N,  etc. 

In  contrast  to  Theorem  3.1  we  cannot  claim  here  that  K 1 when 
h(A)  0.  But  all  experimental  experience  so  far  indicates  that  K never 
exceeds  a value  of  two  for  a large  class  of  problems. 

Theorem  3. 2 holds  no  matter  which  constant  factors  are  used  in  the 
definition  of  the  quantities  (3.22) , (3.24) . But  the  powers  of  h(qj) 
occurring  in  these  quantities  are,  of  course,  essential.  The  particular 
constants  chosen  in  (3.22) , (3.24)  arose  from  an  asymptotic  analysis  of  the 
case  of  a uniform  mesh  A. 

We  illustrate  Theorem  3.2  on  the  following  numerical  example: 

M °~airax-  + ^ = 0.  X € 2 = (0,1)  X (0,1)  c r2 

X 2m  d.X2 

(3.27) 

u = g,  X 6 


The  boundary  condition  is  chosen  to  provide  for  a specific  exact  solution 
Uq.  In  particular,  we  consider  the  following  two  cases: 


1/2 

(I)  a = 0,  Uq  = Re(z-ZQ)  , z - Xj^  + ix2,  Zq  * 1.03, 
(II)  a - 1,  Uq  = Re(w-WQ)^/^,  w * y^  + iy2,  Wq  = yj  + iy® 


where 
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and 


Five  different  meshes  are  used,  denoted  by  (a) , (b) , (c) , (d) , (e) . The 
meshes  (a),  (b),  (c)  are  regular  partitions  of  the  unit  square  with  step- 
size  h ^ J » respectively.  The  other  two  meshes  (d)  and  (e)  are 
shown  in  Figure  4. 


Mesh  (d) 


Figure  4 


Mesh  (e) 


Conforming  bilinear  elements  were  used.  Table  3 gives  the  corresponding 
results.  Here  "active  nodes"  means  the  nunber  of  degrees  of  freedom,  that 
is,  the  size  of  the  corresponding  system  of  linear  equations;  the  other 
notation  is  the  same  as  that  used  in  Table  2 for  the  one -dimensional  case. 

Once  again,  the  a-posteriori  estimates  are  clearly  very  accurate  and  reliable. 
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Table  3 


Case 

active 

nodes 

Jumber  of 
irreg. 
nodes 

elem. 

Energy  < 
abs. 

llellE  ; 

jrror 
in  1 

E 

Effect . 
Const . 

0 

Ratio  of 
error  ind. 

CO 

1(a) 

9 

C 

16 

.0867 

13.51 

.917 

1.17(3) 

1(b) 

49 

0 

64 

.0492 

7.68 

.914 

2.54(4) 

1(c) 

225 

0 

256 

.0253 

3.95 

.921 

3.43(5) 

1(d) 

21 

1 

4 

28 

.0534 

8.35 

.922 

1.31(3) 

1(e) 

33 

8 

40 

.0271 

5.79 

.952 

6.22(3) 

11(a) 

9 

0 

16 

.0677 

11.28 

.683 

7.09(2) 

11  (b) 

49 

0 

64 

.0384 

6.40 

.711 

3.98(3) 

11(c) 

225 

0 

256 

.0201 

3.35 

.736 

1.57(4) 

11(d) 

21 

4 

28 

.0418 

6.97 

.731 

3.03(2) 

11  (e) 

33 

8 

40 

.0294 

4.90 

1 

.805 

2.14(2) 

■»  y-  « '.  •"■■■*' 
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4.  Optimal  Meshes 

4.1  One- Dimensional  Problems 

Once  again  we  consider  the  model  problem  (3.1).  The  theory  of  a- 
posteriori  estimates  sketched  in  Section  3.1  provides  a basis  for  the 
construction  of  optimal  meshes  (3.2)  in  this  case. 

For  this,  we  call  a partition  (3.2)  a (5, m) -partition  if  there  exists 
a continuous,  increasing  function  ^ on  [0,1]  which  satisfies 


(4.1) 


?(Xj)  = ^ , j = 0,l,...,m(A) 


and  is  twice  continuously  differentiable  on  each  subinterval  (x^  2’ 

j = l,...,m.  Any  partition  A is  a (5, m) -partition  for  the  piecewise  linear 

function 


(4.2) 


5(x)=i^  + ^ (x-Xj_j),  x/j  < X < x^,  j = l,...,m 


The  error  corresponding  to  the  (?,m) -partition  A shall  be  denoted 


by 


e(?,m), 


It  now  turns  out  that  the  (Cg,m) -partition  Aq  defined  by 
X 1 

(4.3)  5o(x)  = Yg  / [a(t)u”(t)^]^/^dt,  ^ [a(t)u'(J(t)2]^/^dt 

is  asynptotically  optimal  in  the  following  sense  ([6]): 

Theorem  4.1:  Suppose  that  Ug(x)  has  only  single  roots  in  [0,1], 
then  the  following  statements  hold: 

(a)  The  (Cg,m) -partition  Ag  of  (4.3)  satisfies  the  K-regularity 
condition  (3.4)  with  k = 5/3  for  all  m. 
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(b)  The  error  for  the  partition  Aq  is  given  by 

(4.4)  ||e(Co,m)||g ^ (l-K)(h(Ag)) , as  hCAp)  - 0, 

12m  Yq 

where  the  constant  in  the  bound  of  the  0-term  depends  on  Cq  but  not  on  m. 

(c)  Let  P be  a family  of  5/3-regular  (^,m) -partitions.  Then 


(4.5) 


l|e(4Q,m)|| 


for  any  partition  in  P with  C / ?q. 

Part  (c)  of  this  theorem  shows  that  the  partition  Aq  of  (4.3)  is 
asyn^^totically  optimal  in  the  sense  that  for  any  other  partition  in  P the 
error  is  larger  for  sufficiently  large  m.  Accordingly,  we  call  Aq  the 
asyn5)totically  optimal  mesh. 

The  error  indicators  (A)  of  (3.7)  constitute  a simple  tool  for  the 
construction  of  meshes  with  errors  that  are  asymptotically  equal  to  ||e(^Q,m)||g. 
This  is  the  content  of  the  following  result  (see  again  [ 6 ] ) : 

Theorem  4.2;  Si^pose  that  Uq(x)  has  only  sinqjle  roots  in  [0,1]. 

Let  A be  any  partition  (3.2)  for  which  there  is  a constant  p such  that 


(4.6)  6j  (A)  = p(l+0(h(A)'^)) , K = 1/12,  as  {i(A)  -♦  0 . 

Then  A satisfies  (3.4)  with  k = 5/3  and 


(4.7) 


||e(A)||£  < ||e(AQ)||g(l+0(h(A))''),  as  h(A)  - 0 . 


The  condition  (4.6)  means  that  the  error  indicators  should  be  asympto- 
tically equal.  This  provides  the  basis  for  the  construction  of  partitions 
A satisfying  (4.7).  A natural  approach  is  here  the  adaptive  mesh-refinement 
algorithm  discussed  in  the  next  section. 

It  can  also  be  shown  (see  [6])  that  the  error  for  the  (?Q,m) -partition 
is  very  stable.  Nbre  precisely,  if  A is  any  (C,m) -partition  with 
5 + Xq  then 

(4.8)  l|e(C,m)||g  = l|e(^Q,m)||g(l+0(x2))  as  X ^ 0 . 

On  the  other  hand,  the  distribution  of  the  nodal  points  of  the  optimal 
partition  is  not  particularly  stable.  In  other  words,  we  should  concentrate 
on  making  the  error  indicators  nearly  equal  and  not  concern  ourselves 
with  an  accurate  computation  of  Aq  itself. 

Numerical  results  illustrating  these  results  were  already  included  in 
Tables  1 and  2 of  Section  3.1.  In  particular,  in  the  second  part  of  Table  2 
the  00-values  are  nearing  one  and  as  expected  the  values  of  E and  Eq  are 
very  close. 

4.2  Higher-Dimensional  Problems 

For  problems  in  two-  or  higher  dimensions  there  exists  as  yet  no 
theory  of  the  extent  sketched  in  Section  4.1  for  the  one -dimensional  case. 

In  [ 3 ] it  is  shown  that  here  again  asynptotic  optimality  of  the  mesh  calls 
for  the  (asyjujtotic)  equality  of  the  error  indicators.  This  provides  the 
basis  for  the  following  type  of  adaptive  mesh- refinement  algorithm. 
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Suppose  that  the  error  indicators  have  an  asymptotic  behavior 
of  the  form 

(4.9)  ^i  “ ^i*  as  h^  -»■  0 . 

If  any  element  with  corresponding  value  was  generated  by  sub- 
dividing an  element  in  a previous  mesh  with  value  then  (4.9) 

suggests  that  the  (worst)  indicator  after  dividing  will  be  approximately 

(4.10)  = (e.)2/e?^‘^  . 

Practical  experience  has  shown  that,  in  general,  this  prediction  is  very 
satisfactory. 

Clearly,  we  should  refine  only  those  elements  which  have  an  indicator 
above  the  largest  predicted  new  indicator  in  the  new  mesh.  In  order  to  , 
start  this  process,  all  elements  in  the  basic  mesh  are  refined  in  the  first 
step.  Hence  we  have  the  following  refinement  algorithm: 

1.  cut  :=  0 

2.  If  "current  mesh  L is  the  basic  mesh"  then  go  to  4 

3.  For  "each  element  t in  A"  do 

3.1  Confute  error  indicator  e of  element  x 

, new  _ 2 / old 

o*  it  6 “6/6 

3.3  If  > cut  then  cut  :■  e”®^ 

4.  For  "each  element  x in  A"  do 

4.1  If  6 > cut  then  "subdivide  x and  for  each  new 
element  set  :=  e 

Implementation  details  will  be  discussed  elsewhere. 
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5.  Future  Trends 

5.1  Computational  Goals 

In  the  development  of  typical  mathematical  library  programs,  the  aim 
is  to  achieve- -whenever  possible--a  prescribed  and  typically  high  accuracy. 
Here  "accuracy"  refers  usually  to  the  size  of  the  errors  introduced  by  the 
particular  numerical  method,  and  the  program  is  reliable  if  its  overall 
failure  rate  is  low. 

For  the  general-purpose,  finite -element  software  used,  say,  in  struc- 
tural engineering  the  requirements  are  rather  different.  The  user  needs  to 
obtain  results  which  predict,  with  an  acceptable  degree  of  reliability  and 
accuracy,  the  actual  behavior  of  the  mechanical  structure  at  hand.  In  other 
words,  the  terms  "accuracy"  and  "reliability"  are  used  in  a much  broader 
context  than  before.  It  is  clearly  inadequate  to  assess  the  results  solely 
on  the  basis  of  the  errors  introduced  by  the  numerical  calculations.  For 
example,  an  earlier  decision  to  use  a plate  model  rather  than  a full  three- 
dimensional  formulation  may  result  in  final  stress  values  that  differ  by 
10-20%  from  each  other.  In  contrast  to  this,  the  numerical  errors  are 
entirely  negligible.  In  other  words,  the  accuracy  of  the  numerical  results 
should  be  a measure  of  their  deviation  from  the  solution  of  a "higher" 
matliematical  model.  Moreover,  these  results  can  only  be  considered  reliable 
if  changes  in  the  reference  model  and,  more  generally,  in  the  entire  sequence 
of  steps  from  the  problem  formulation  to  the  final  results  have  relatively 
small  effects. 
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Generally,  the  confutations  involved  in  finite-element  applications  are 
ver)'  extensive  and  highly  demanding  on  the  available  computer  resources,  not 
to  mention  the  corresponding  demands  on  the  user  personnel.  In  most  cases 
the  user  is  forced  to  set  a maximum  allowable  cost  for  any  one  computational 
analysis.  Ideally  then  the  goal  is  to  achieve  optimal,  acceptable  accuracy 
(in  the  above  sense)  within  the  prescribed  cost  range.  Here  it  is  essential 
that  often  relatively  low  accuracies --often  in  the  range  of  5-10l--are 
entirely  satisfactory. 

Without  doubt  the  design  goal  of  producing  reliable  results  of  optimal 
accuracy  within  some  prescribed  cost  range  represents  a serious  challenge 
and  is  not  realized  in  today's  finite- element  software.  However,  our  own 
exi^erience  has  shown  that  this  is  not  an  unrealistic  goal  (see,  e.g. , [3  ], 

[4  ]>  [13] » arid  it  appears  also  that  it  corresponds  well  with  the  trends 
that  are  now  beginning  to  emerge  in  the  design  of  the  next  generation  of 
such  software. 

The  principal  tool  for  the  realization  of  this  design  goal  is  the 
availability  of  relatively  easily  confutable  error  estimates.  On  some  simple 
model  problems  we  have  illustrated  in  the  previous  two  sections  that  a- 
posteriori  error  estimates  are  feasible  and  not  very  costly  to  compute. 
Moreover,  these  estimates  are  very  reliable.  The  approaches  presented 
here  extend  readily  to  much  more  general  cases  including,  for  instance, 
the  typical  elasticity  problems  in  structural  analysis 

and  even  nonlinear  problems.  Moreover,  the  estimates  can  also  be  developed 
for  norms  other  than  the  energy  norm.  Of  course,  the  formulas  for  the  error 
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indicators  are  then  different,  but  the  structure  of  the  estimates  remains 
the  same. 

I 

Tlie  a-posteriori  estimates  are  developed  on  the  basis  of  a general 
definition  of  the  finite  element  method  in  terms  of  a weak  mathematical 
formulation  involving  bivariate  forms  on  suitable  function  spaces  (see, 
e.g.,  [13]).  This  suggests  that  finite  element  software  should  no  longer 
be  designed  for  specific  classes  of  applications  but  instead  should  consti- 
tute more  general  "finite  element  solvers"  based  on  entire  classes  of  such 
weak  formulations.  This  will  provide  the  needed  closer  interaction  with 
the  theoretical  foundations  and,  at  the  same  time,  it  will  open  up  a much 
wider  applicability  of  the  software. 

In  order  to  simplify  the  use  of  the  software,  the  user  should  not 
have  to  understand  the  internal  operation  or  to  make  difficult  a-priori 
decisions  about  the  desired  computations.  For  this  the  user  interfaces 
should  be  strictly  separated  into  flexible  pre-  and  post -processors  which 
adapt  the  systan  to  various  classes  of  applications. 

The  need  for  reducing  the  range  of  decisions  asked  of  the  user  in 
today's  finite  element  software  requires  the  introduction  of  extensive 
adaptive  techniques  into  the  computational  process.  These  techniques  are 
also  needed  to  meet  the  goal  of  achieving  an  optimal  solution  within  a 
prescribed  cost  range.  The  adaptive  control  of  the  computation  in  turn 
requires  the  availability  of  reliable  error  estimates  at  all  stages  of  the 
calculation.  In  other  words,  here  is  a further  reason  for  the  introduction  ^ 

of  the  a-posteriori  estimates  discussed  in  the  previous  sections.  1 

I 

I 
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In  connection  with  the  optimality  requirement  we  have  to  distinguish 
between  the  design  of  a (nearly)  optimal  mesh  and  the  achievement  of  a 
(nearly)  optimal  a-posteriori  error  estimate  for  the  solution.  As  we  saw 
before,  the  optimal  error  (for  a given  degree  of  freedom)  is  rather  stable 
while  this  is  not  true  for  the  optimal  mesh.  Thus  we  should  certainly 
concentrate  on  achieving  a (nearly)  optimal  error  within  cur  given  cost 
range  and  not  overplay  the  search  for  the  optimal  mesh.  The  discussions  of 
Section  4 show  that  we  have  a rather  siirple  criterion  for  the  (near)  opti- 
mality of  the  error  in  the  near-equality  of  the  error  indicators.  It 
provides  the  basis  of  single  and  yet  highly  effective  adaptive  mesh-refine- 
ment  algorithms  (see,  e.g.,  [3],  [13]).  Since  these  algorithms  aim  to  keep  all 
error  indicators  as  close  together  as  possible,  all  resulting  meshes--after 
some  initial  phase- -have  a fairly  optimal  error  for  the  degrees  of  freedom 
they  incorporate.  Thus  we  may  stop  the  process  whenever  either  a prescribed 
error  tolerance  lias  been  achieved  or  the  given  computational  cost  range  has 
been  exceeded  (see  [3], [4]). 

5. 2 A Prototype  Adaptive  Finite  Element  Solver 

As  stated  before,  the  general  design  goals  for  future  finite  element 
software  still  represent  serious  challenges  and  are  not  yet  incorporated 
into  any  operational  programs.  However,  as  we  said,  they  are  within  the 
realm  of  practicality. 

In  substantiation  of  this  claim  an  experimental  finite  element  system 
is  being  designed  by  several  of  us  at  the  University  of  Maryland  that  meets 
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the  following  four  design  criteria: 

(a)  The  system  constitutes  an  applications- independent  finite  element 
solver  for  a certain  class  of  two-dimensional,  linear,  elliptic 
problans  based  on  a weak  mathematical  formulation; 

(b)  it  incorporates  extensive  adaptive  approaches  to  minimize  the 
critical  decisions  demanded  of  the  user; 

(c)  it  uses  a-posteriori  error  estimates  of  the  type  discussed  above 
to  control  the  adaptive  approach  and  to  provide  a solution  with 
a (near)  optimal  error  within  a prescribed  cost  range; 

(d)  it  takes  advantage  of  modem  advances  in  the  design  of  data 
structures  and  software  systems,  and,  in  particular,  it  is  designed 
to  use  the  natural  parallelism  and  modularity  of  the  finite  element 
method  to  increase  the  size  of  the  practically  solvable  problems. 

The  general  design  of  the  system  has  been  described  in  [13] . It  is 
being  iii^Dlemented  not  as  a production  system  but  as  an  experimental  system 
for  the  evaluation  of  the  various  new  ideas  in  it.  Accordingly,  extensive 
provisions  for  evaluating  the  performance  are  incorporated  in  it. 

Beyond  the  new  approaches  reflected  in  the  criteria  (a) , (b)  and  (c) 
already  discussed  in  the  previous  section,  the  introduction  of  parallelism 
appears  to  be  a particularly  novel  aspect  of  the  design.  This  parallelism 
is  on  the  procedural  level  rather  than  the  instructional  level  because  there 
the  expected  payoff  is  much  greater.  Thus,  parallelism  is  specified  in  terms 
of  processes  which  are  autonomous  units  with  their  own  programs  and  data. 
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The  processes  run  in  parallel  and  communicate  asynchronously  in  a limited 
and  highly  structured  manner. 

The  natural  parallel  process  structure  for  the  finite  element  system 
derives  from  the  familiar  substructure  analysis  in  engineering  design. 

In  principle  the  substructure  segmentation  of  the  data  and  processing  is 
the  same  as  that  used  in  large  industrial  applications  except  that  ours 
is  essentially  automatic  and  more  flexible,  while  current  finite  element 
systems  denond  that  their  users  create  a rigid  segmentation  at  the  manual 
or  even  managerial  level. 

In  the  design  the  local  data  of  a process  contain  almost  all  the  in- 
formation needed  for  the  computations  of  that  process.  Thus  the  process 
structure  induces  on  the  finite  element  data  a segmentation  which  is  natural 
to  the  problem  and  provides  a basis  for  intelligent  storage  management  in 
the  environment  of  a single  large  computer.  At  the  same  time  the  use  of 
parallel  processes  makes  it  possible  to  apply  multiple  processors  effec- 
tively- -hopefully  for  significant  gains  in  speed. 
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